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[j, • The Navier-Stokes transport coefficients for binary mixtures of smootli inelastic liard disks or 

splieres under gravity are determined from tiie Boltzmann kinetic tiieory by application of the 
Chapman-Enskog method for states near the local homogeneous cooling state. It is shown that the 
Navier-Stokes transport coefficients are not affected by the presence of gravity. As in the elastic 
case, the transport coefficients of the mixture verify a set of coupled linear integral equations that 
ri , are approximately solved by using the leading terms in a Sonine polynomial expansion. The results 

^ ■ reported here extend previous calculations [V. Garzo and J. W. Dufty, Phys. Fluids 14:1476-1490 

D ' (2002)] to an arbitrary number of dimensions and provide explicit expressions for the seven Navier- 

Stokes transport coefficients in terms of the coefficients of restitution and the masses, composition, 
and sizes of the constituents of the mixture. In addition, to check the accuracy of our theory, the 
inelastic Boltzmann equation is also numerically solved by means of the direct simulation Monte 
Carlo method to evaluate the diffusion and shear viscosity coefficients for hard disks. The compar- 
^ , ison shows a good agreement over a wide range of values of the coefficients of restitution and the 

parameters of the mixture (masses and sizes) . 
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I. INTRODUCTION 



f^ ' The simplest model for a granular fluid is a system composed by smooth hard spheres or disks with inelastic 

^ . collisions. The only difference from the corresponding model for normal fluids is the loss of energy in each binary 
fS| collision, characterized by a (constant) coefflcient of normal restitution. For a low density gas, the Boltzmann kinetic 

VO ! equation conveniently modified to account for inelastic collisionsii^i^ has been used in recent years as the starting point 
*^ ■ to derive the hydrodynamic-like equations of the system. Thus, assuming the existence of a normal (hydrodynamic) 
^r ' solution for sufficiently long space and time scales, the Chapman-Enskog (CE) method^ has been applied to solve the 
^^i Boltzmann equation to Navier-Stokes (NS) order and get explicit expressions for the transport coefficients. While 
f^ this goal has been widely covered in the case of a monocomponent gaSf^ much less is known for systems composed by 
"T-*, , grains of different masses, sizes, and concentrations (granular mixtures). 

C^ ' Needless to say, the determination of the NS transport coefficients for a multicomponent granular fluid is much 

U , more complicated than for a single granular system. Many attempts to derive these coefficients® have been carried 

I ■ out by means of the CE expansion around Maxwellians at the same temperature T for each species. The use of this 

■ O ' distribution can only be considered as acceptable for nearly elastic systems where the assumption of the equipartition 

^ . of energy still holds. In addition, according to this level of approximation, the inelasticity is only accounted for by 

the presence of a sink term in the energy balance equation and so the expressions of the NS transport coefficients are 

the same as those obtained for elastic collisions. However, as the dissipation increases, different species of a granular 

mixture have different partial temperatures Ti and consequently, the energy equipartition is seriously broken (Ti ^ T) . 

The failure of energy equipartition in granular fluids^'* has also been confirmed by computer simulations- and even 

5J] ' observed in real experimentsiS, of agitated mixtures. All the results show that deviations from equipartition depend 

Ci I on the mechanical differences between the particles of each species and the coefficients of restitution of the system. 

Given that the inclusion of nonequipartition effects increases the level of complexity of the problem, it is interesting 

from a practical point of view to assess the influence of this effect on the transport properties of the system. If the 

NS transport coefficients turned out to be quite sensitive to nonequipartition, the predictions made from previous 
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theories^ should be reexamined by theories that take into account the nonequipartition of energy. For this reason, 
although the possibility of nonequipartition was already pointed out many years agOfii a careful study of its influence 
on transport has only been carried out recently. In this context, Garzo and Duftj>i^ have developed a kinetic theory for 
a binary granular mixture of inelastic hard spheres at low density which accounts for nonequipartition effects. Their 
results show that in general the consequences of the temperature differences for the transport coefficients are quite 
significant, especially for strong dissipatioufi^^ It is important to remark that the expressions derived in Ref. Il2 for 
the Navier-Stokes transport coefficients do not limit their application to weak inelasticity. In fact, the results reported 
in this paper include a domain of both weak and strong inelasticity, 0.5 < a < 1, where a is the (common) coefficient 
of restitution considered. The accuracy of these theoretical predictions (based on a Sonine polynomial expansion) 
has been confirmed by Monte Carlo simulations of the inelastic Boltzmann equation in the cases of the diffusion 
coefficient^ and the shear viscosity coefficient of a mixture heated by an external thermostat J^ Exceptions to this 
good agreement are extreme mass or size ratios and strong dissipation, although these discrepancies are mitigated 
in part if one retains more terms in the Sonine polynomial expansionJ^ For small dissipation, the results derived 
by Garzo and Duftj»i^ agree with those recently obtained by Serero et alr^ in the first order of the order parameter 

The CE method solves the Boltzmann equation by expanding the distribution function of each species /^(r, v, f) 
around the local homogeneous cooling state (HCS).^ This state plays the same role for granular gases as the local 
equilibrium distribution for a gas with elastic collisions. Given that the form of the distribution function f!> of the 
HCS is not exactly known, one usually considers the first correction to a Maxwellian at the temperature for that 
species, namely, a polynomial in velocity of degree four (leading Sonine correction). However, the results derived 
for hard spheres clearly show that the influence of these non-Gaussian contributions to the transport coefficients 
are in general negligible, except in the case of the heat flux for quite large values of dissipation.^** Accordingly, a 
theory incorporating the contributions coming from the deviations of the HCS from its Gaussian form does not seem 
necessary in practice for computing the NS transport coefficients of the mixture. 

The objective of this paper is twofold. First, given that the results reported in Ref. |l^ are limited to hard spheres, 
we extend here this derivation to an arbitrary number of dimensions d. This goal is not only academic since, from 
a practical standpoint, many of the experiments reported for flowing granular materials have created (quasi) two- 
dimensional systems by confining grains between vertical or on a horizontal or tilted surface, enabling data collection 
by high-speed video. *^ Regarding computer simulations, most of them consider hard disks to save computer time and 
memory. For these reasons, it would be desirable to provide experimentalists and simulators with theoretical tools 
to work when studying problems both in two and three dimensions. In addition, apart from its practical interest, 
it is also interesting from a fundamental view to explore what is the influence of dimensionality on the dependence 
of the transport coefficients on dissipation. As a second target, we want also to present a simplified theory with 
explicit expressions for the transport coefficients. As the algebra involved in the calculations of Ref. [13 is complex, 
the constitutive relations for the fluxes were not explicitly displayed in this paper. Although the work carried out here 
involves complex algebra as well, the use of Maxwellians at different temperatures for the distribution functions of 
each species in the reference state allows us to explicitly obtain expressions for the seven relevant transport coefficients 
of the mixture in terms of the mechanical parameters of the system: masses, sizes, composition and coefficients of 
restitution. To assess the degree of accuracy of our (approximated) expressions, we have also performed Monte Carlo 
simulations for the diffusion and the shear viscosity coefficients for hard disks {d = 2). As shown below, the good 
agreement found between the results derived in this paper with computer simulations justifies this simplification and 
allows one to obtain more simplified forms of the transport coefficients. 

The plan of the paper is as follows. In Sec.|TTl the inelastic Boltzmann equation and the corresponding hydrodynamic 
equations are recalled. The CE expansion adapted to the inelastic binary mixtures is formulated in Sec. lIIII Assuming 
that gradients and dissipation are independent parameters, the Boltzmann equation is solved by means of an expansion 
in powers of the spatial gradients around the local HCS distribution /^ . It is shown that the use of the local HCS as 
the reference state is not an assumption of the CE method but a consequence of the exact solution to the Boltzmann 
equation in the zcroth-order approximation. Section llVl deals with the expressions for the NS transport coefficients. As 
in the case of elastic collisions, these coefficients are the solutions of a set of coupled linear integral equations which 
involves the (unknown) distributions /^ . The integral equations are solved by considering two approximations: 
First, /j is replaced by its Maxwellian form at the temperature Ti, and second, only the leading terms in a Sonine 

polynomial expansion of the first-order distribution /^ are retained. Technical details of the calculations carried out 
here are given in Appendices El [B] and [C] A comparison with previous results'^- based on the use of Maxwellians 
at the same temperature T as a ground state is also illustrated in Sec. II V[ showing significant discrepancies between 
both descriptions at moderate dissipation. Section |V] is devoted to the numerical solutions of the Boltzmann equation 
by using the direct simulation Monte Carlo (DSMC) methodic in the cases of the diffusion D and shear viscosity ?/ 
coefficients for hard disks. To the best of our knowledge, this is the first time that the NS shear viscosity of a granular 



binary mixture at low density has been numerically obtained from the DSMC method. The paper is closed in Sec. IVII 
with a brief discussion of the results presented in this paper. 

II. BOLTZMANN EQUATION AND CONSERVATION LAWS 

Consider a binary mixture composed by smooth inelastic disks {d = 2) or spheres (d = 3) of masses mi and m2, 
and diameters cti and (T2 . The inelasticity of collisions among all pairs is characterized by three independent constant 
coefficients of restitution an, a22, and ai2 ~ a2i, where a^ < 1 is the coefRcient of restitution for collisions between 
particles of species i and j. The mixture is in presence of the gravitational field so that each particle feels the action 
of the force F.^ = ni^g, where g is the gravity acceleration. In the low density regime, the distribution functions 
/i(r,v;t) (z = 1,2) for the two species are determined from the set of nonlinear Boltzmann equations^ 

at+vV + g-— j/,(r,v,i)=^Jy[v|/,(i),/,W] , (1) 

where the Boltzmann collision operator Jij [v|/i, fj] is 

J^3 b^ilftJj] = crff^ / dv2 / daQia ■ gi2)(S^ • gi2) 

x[ar.2/,(r,vl,t)/,(r,v^,i)-/,(r,vi,t)/,(r,V2,i)] . (2) 

In Eq. ([2), d is the dimensionality of the system, aij — (at + Uj) /2, a is an unit vector along the line of centers, Q is 
the Heaviside step function, and gi2 = vi — V2 is the relative velocity. The primes on the velocities denote the initial 
values {v'j^,V2} that lead to {vi,V2} following a binary (restituting) collision: 

v'l = vi - ^j, (1 + a-^) (ct • gi2)CT, 

V2 = V2 + ^liJ (1 + fty^) (S • gl2)ff, (3) 

where /Xy = mi/ {m,i + m,j). The relevant hydrodynamic fields are the number densities n^, the flow velocity u, and 
the temperature T . They are defined in terms of moments of the distributions fi as 

2 
"•j = / dvf,{v) , pu = ^m^ / dvv/i(v) , (4) 

2 2 

nT = p = ^n,T,-^^ /dvl/V.(v), (5) 

where V = v — u is the peculiar velocity, n — ni + n2 is the total number density, p — mini + m2n2 is the total mass 
density, and p is the pressure. Furthermore, the third equality of Eq. ([5|) defines the kinetic temperatures Ti for each 
species, which measure their mean kinetic energies. 

The collision operators conserve the particle number of each species and the total momentum but the total energy 
is not conserved: 

I dvJ.,,[v\f,J,] ^ , (6) 

22 „ 

^ ^ m, / dvvJ,, [v|/„ /,] = , (7) 

i=l 3 = 1 ■' 

2 2 „ 

Y.Y.'^-j d^y^J^o [v|/^, /,] = -dnTC , (8) 

i=l j = l 



where ( is identified as the total "cooling rate" due to inelastic collisions among all species. At a kinetic level, it is 
also convenient to introduce the "cooling rates" Q for the partial temperatures T^. They are defined as 



2 2 



C» = EC^^ = -i^.ilfd^y'M^\hJA^ (9) 

where the second equality defines the quantities Qj . The total cooling rate C can be written in terms of the partial 
cooling rates C,i as 

2 

C = r-i^x,T,C., (10) 

where Xi = nijn is the mole fraction of species i. 

From Eqs. ([4])-([8|), the macroscopic balance equations for the mixture can be obtained. They are given by 

An. + n,V-u+^^^ =0, (11) 

Au + p-^V-P^g, (12) 

Ar--y ^^^ + :^(v-q+P: Vu) = -cr. (13) 

n ^-^ TJii an 

i—l 

In the above equations, I?t = 9t + u • V is the material derivative, 



j. =m,ydvV/,;(v) (14) 



is the mass flux for species i relative to the local flow, 

P = E™' / dvVV/,(v) (15) 

i=i •' 

is the total pressure tensor, and 

ci = Y^'^ I d^V'Y n{.) (16) 

i=l "' 

is the total heat flux. 

The macroscopic balance equations (|lip -(fT3 |) are not entirely expressed in terms of the hydrodynamic fields, due 
to the presence of the cooling rate C, the mass flux j^, the heat flux q, and the pressure tensor P which are given as 
functionals of the distributions fi. However, it these distributions can be expressed as functionals of the hydrodynamic 
fields, then the cooling rate and the fluxes also will become functional of the hydrodynamic fields through Eqs. Q 
and (|14p - (|16p . Such expressions are called constitutive relations and they provide a link between the exact balance 
equations and a closed set of equations for the hydrodynamic fields. This hydrodynamic description can be derived 
by looking for a normal solution to the Boltzmann kinetic equation. A normal solution is one whose all space and 
time dependence of the distribution function fi occurs through a functional dependence on the hydrodynamic fields, 

/,(r,v,t) = /,[v|a:i(r,t),p(r,i),r(r,i),u(r,i)] ■ (17) 

As in previous worksji^iii^ we have taken the set {a;i,p, T, u} as the d + 3 independent fields of the two-component 
mixture. These are the most accessible fields from an experimental point of view. The determination of this normal 
solution from the Boltzmann equation ([T]) is a very difficult task in general, unless the spatial gradients are small. In 
this case, the CE method gives an approximate solution. 



III. CHAPMAN-ENSKOG SOLUTION 

The CE method is a procedure to construct an approximate normal solution. It is perturbative, using the spatial 
gradients as the small expansion parameter. More specifically, one assumes that the spatial variations of the hydrody- 
namic fields rii, u, p, and T are small on the scale of the mean free path. For ordinary gases this can be controlled by 
the initial or boundary conditions. It is more complicated for granular gases, since in some cases (e.g., steady states 
such as the simple shear flow problem^S.) the boundary conditions imply a relationship between dissipation and some 
hydrodynamic gradient. As a consequence, there are examples for which the NS approximation is restricted to the 
quasi-elastic limiti^ Here, we also assume that the spatial gradients are independent of the coefficients of restitution 
so that, the corresponding NS order hydrodynamic equations apply for small gradients but they are not limited a 
priori to weak inelasticity. It must be emphasized that our perturbation scheme differs from the one recently carried 
out by Serero et alt^ where the CE solution is given in powers of both the hydrodynamic gradients (or equivalently, 
the Knudsen number) and the degree of dissipation e^ = 1 — af . The results provided in Ref. [iTJ only agree with 
our results in the quasielastic domain (small 6^). Moreover, in the presence of an external force it is necessary to 
characterize the magnitude of the force relative to gradients as well. As in the elastic casej^ it is assumed here that 
the magnitude of the gravity field is at least of first order in perturbation expansion. 

For small spatial variations, the functional dependence (J17p can be made local in space through an expansion in 
gradients of the hydrodynamic fields. To generate it, fi is written as a series expansion in a formal parameter 5 
measuring the nonuniformity of the system. 



_ Ao) 



fr = fr+sfr+s'fr + 



(1) 



f(2) 



(18) 



where each factor of 5 means an implicit gradient of a hydrodynamic field. The local reference states /^ are chosen 
such that they verify Eqs. ^ and ((5|), or equivalently, the remainder of the expansion must obey the orthogonality 
conditions 



dv 



/.(v)-/f (v) 



= 



(19) 
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1=1 



dvv 



/.(v)-/f(v) 



(20) 



E rrij 
~9" 



dvV^ 



/.(v)-/f"^(v) 



= 



(21) 



The time derivatives of the fields are also expanded as dt = 9, 



(0) 



eo: 



(1) 



The action of the operators d^ 



(k) 



can be obtained from the balance equations (|lip ~ (|13p when one takes into account the corresponding expansions for 
the fluxes and the cooling rate. This is the usual CE method"' for solving kinetic equations. The main difference in 
the case of inelastic collisions is that the reference state has a time dependence associated with the cooling that is 
not proportional to the gradients. As a consequence, terms from the time derivative d^ are not zero. In addition, 

(k) 

the different approximations /^ are well-defined functions of the coefficients of restitution a^, regardless of the 
applicability of the corresponding hydrodynamic equations truncated at that order. 



A. Zeroth-order approximation 

To zeroth order in the gradients, Eq. ([T|) becomes 



(22) 



where use has been made of the fact that gravity is assumed to be of first order in the uniformity parameter S. The 
balance equations to this order give 



af ^xi = di"\i = 0, T-'d^°h = p-iaf V = -C<"^ 



(23) 



where C^°^ is determined by Eqs. ([9]) and (fTO|) to zeroth order in the gradients. Since f> is a normal solution, then 
the time derivative in Eq. (j22p can be written as 



5(0)^(0) ^ _^(o) (^Tdr+pd,) /(°) = ic(°^^ • (v/f ^ 



(24) 



The second equality in Eq. (j24p follows from dimensional analysis which requires that the dependence of /^ on p 
and T is of the form 



f^'\v)^x,^v,'<^dV/vo), 



(25) 



where vo{t) = y/2T{mi + m2)/'inim2 is a thermal velocity defined in terms of the temperature T(t) of the mixture 

and $i is a dimensionless function of the reduced velocity V/vq. The dependence of /| ^ on the magnitude of V 
follows from the isotropy of the zeroth-order distribution with respect to the peculiar velocity. Thus, the Boltzmann 
equation at this order reads 






(26) 



Equation ([26]) has the same form as the Boltzman equation for a strictly homogeneous state. The latter is called 
the homogeneous cooling state (HCS).^ Here, however, the state is not homogeneous because of the requirements 
p9)) -- (|2T|) . Instead it is a local HCS. It must be emphasized that the presence of this local HCS as the ground or 
reference state is not an assumption of the CE expansion but rather a consequence of the kinetic equations at zeroth 
order in the gradient expansion. 

The local HCS distribution is the solution of the Boltzmann equation (pS)) . However, its explicit form is not exactly 
known even in the one-component case.^^ An accurate approximation for the zeroth-order solution f^ can be obtained 

by using low order truncation of a Sonine polynomial expansion. The results show that in general, f^ is close to 
a Maxwellian at the temperature for that species. Further details of this solution for hard spheres {d = 3) can be 
found in Ref. |7|. An important consequence is that the kinetic temperatures of each species are different for inelastic 
collisions and, consequently the total energy is not equally distributed between both species (breakdown of energy 
equipartition). This violation of energy equipartition has been confirmed by computer simulation studies^ as well as 

by real experiments)^ The condition that /j is normal in the sense of Eq. (jl7[) (namely, it depends on time only 
through its functional dependence on T and p) implies that the ratio Ti/T = ^i{xi) depends on the hydrodynamic 
state through the concentration xi. 

The dependence of the temperature ratio 7 = 71/72 = T1/T2 on the parameters of the mixture is obtained by 
requiring that the partial cooling rates Q must be equal,^ i.e.. 



cf) = cf = c(°^ 



(27) 



(0) 



These partial cooling rates are nonlinear functionals of the distributions /^ , which are not exactly known. However 
to get the temperature ratio, they can be well estimated by using Maxwellians at different temperatures: 



/f^(V)^/,M(V) 



/ TOi 



V27rT, 



d/2 



exp 



vriiV 



(28) 



In this approximation, one get; 



,22 



j=l V2/ j=l 

x(l + a,,)[l-^(l + a,,)^ 



1/2 



(29) 



where 



7* 



E "1, '• 



j=i 



(30) 



It must be remarked that the fact that Ti(i) ^ T2{t) does not mean that there are additional hydrodynamic degrees 
of freedom since the partial temperatures Ti can be expressed in terms of the granular temperature T as 



Tiit) 



7 



l + a;i(7-l) 



nt), T2{t) = 



1 



l + a;i(7-l) 



nt). 



(31) 



Note that the reference Maxwellians ([^5]) for the two species can be quite different due to the temperature differences. 

This contrasts with the ground state considered in more standard derivation a^'^^ where /^ is replaced by a Maxwellian 
defined at the same temperature T, i.e., 



As will show later, the approaches 



/r(v) 

and 



/ rui \d/i 



2T 



(32) 



lead to different results for the NS transport coefficients. 




FIG. 1: Temperature ratio 7 = T1/T2 versus the (common) coefficient of restitution a for an equimolar mixture [xi — 1/2) of 
hard disks {d — 2) with u! = fJ.^^^- Two different vafues of the mass ratio are considered; /i = 4 and fi = 8. The symbols refer 
to DSMC results while the lines represent the theoretical results obtained from the condition (|27|) . 



The solution to Eq. (P7|) gives the temperature ratio T1/T2 for any dimension d as a function of the mole fraction 
xi, the mass ratio /i = 7711/7712, the size ratio uj = cti/(T2, and the coefficients of restitution a^. To illustrate the 
violation of equipartition theorem, in Fig. [1] we plot the temperature ratio versus the coefficient of restitution a for 
hard disks (d = 2) in the case of an equimolar mixture Xi = ^ for two different mixtures composed by particles 
of the same material: ^ = 4, w = 2, and /i = 8, a; = v8- For the sake of simplicity, we have taken a common 
coefficient of restitution a = an — 0,22 — cki2- We also include the simulation data obtained by solving numerically 
the Boltzmann equation by means of the DSMC methodi^ The excellent agreement between theory and simulation 
shows the accuracy of the estimate (1^^ to compute the temperature ratio from the equality of cooling rates (|27p . 
We also observe that the deviations from the energy equipartition increase as the mechanical differences between the 
particles of each species increase. 



IV. NAVIER-STOKES TRANSPORT COEFFICIENTS 



The CE procedure allows one to determine the NS transport coefficients of the mixture in the first order of the 
expansion. The analysis to first order in gradients is similar to the one worked out in Ref. [12 for d = 3. Here, we 
only present the final results with some technical details being given in Appendix[Xl The mass, momentum, and heat 
ffuxes are given, respectively, by 



:(1) 
Jl 



7771777,277 



D^xy- -DpVp 



T 



,(1) 
J2 



(33) 



P, 



(1) _ 



kf. 



pSki-il [ V^Ufc + VfcTt^ 



-jSki"^ ■ u 



(34) 



q(i) ^ -T^D"Vxi - LVp - AVT. 
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(35) 



The transport coefScients in these equations are the diffusion coefficient D, the pressure diffusion coefficient Dp, the 
thermal diffusion coefficient D', the shear viscosity 77, the Dufour coefficient D" , the pressure energy coefficient L, 
and the thermal conductivity A. These coefficients are defined as 



D = 



dm2n 



dvV-Ai, 



(36) 



Dp^-^ I d^Y-B,, 



^''-^Z— ■■ 



(37) 
(38) 



rj = 



id- 



2 

l)(d + 2)^ 'J 



(39) 



i— 1 



(40) 



1 V^ m^ 



^-iEv^-^^v-^" 






(41) 



A-jE^/'ivl'^V.C, 



d^^ 2 



(42) 



As for ordinary gases,"* the unknowns ^(V), Bi(V), Ci(V), and X>i(V) are the solutions of the following set of 
coupled linear integral equations: 



-C(°)(TaT+pap) + A 



Ai + M1A2 = Ai + (^^] [pBi + TCi 



\dx, J 



P,T 



(43a) 



(^''HTdT + pdp) + C2 



A2 + M2A1 = A2 + (^) {pB2 + TC2) 

\ 0^1 Jp,T 



(43b) 



-C^''HTdT+pdp)+Ci-2C 



(0) 



y/-(0) 

Bi+MiB2 = Bi + -^ei, 



(44a) 



<(°) (T5T+pap)+ £2 -2C("^ 



B2+X2Bl=B2 + ^— C2, 



(44b) 



-C(°)(raT+pap) + /:i-ic^°) 






(45a) 



-C(°)(raT + pap) + /:2-2C^"^ 



C2+M2Ci=C2-^B2, 



(45b) 



C^^HtOt + pdp) + Ci 



■D1+M1V2 = Di, 



(46a) 



C^^HtOt + pdp) + C2 



V2+M2'Di ^ D2 



(46b) 



In the above equations, the quantities A^, B^, C^, and D^ are given by Eqs. (JA8|) - (|A11[) . respectively. They depend 
on the local HCS distribution /^ . In addition, we have introduced the linearized Boltzmann collision operators 



CiX = -(^Jn [/f ^ , A] + Jn [A, /(°)] + J,2 [A, /f )] 



(47) 



MiX = -J,2[.fr.X] 



(0) 



(48) 



The corresponding expressions for the operators £2 and AI2 can be easily obtained from Eqs. (|47|) and (|48|) by just 
making the changes 1 ^^ 2. Note that in Eq. (|15)l the cooling rate C^"-* depends on xi explicitly and through its 
dependence on "f{xi). This dependence gives rise to significant new contributions to the integral equations for the 
transport coefficients. Furthermore, the external field does not occur in Eqs. (|43p - (l46l) . This is because the particular 
form of the gravitational force. 



A. Sonine polynomial approximation 



So far, all the results are exact. However, explicit expressions for the NS transport coefficients requires to solve 
Eqs. (|43p - (P5|) as well as the integral equations (pS)) for the reference distributions f^ '. As said before, the results 
obtained in the HCSZ have shown that /} is well represented by its Maxwellian form (P5|) in the region of thermal 
velocities. For this reason and to provide simple but accurate expressions for the transport coefficients, non-Gaussian 
corrections to /,| ' will be neglected in our theory. The full expressions for the transport coefficients in the case 
d = 3 (including non-Gaussian corrections) can be found in Ref. Il2l . It must remarked that while the effect of these 
non-Gaussian corrections on the transport coefficients is not important in the case of the mass flux and the pressure 
tensor, the same does not happen for the heat flux, where the influence of them is not negligible at high inelasticity.— 
With respect to the functions {jA.i,Bi,Ci,'Di), we will expand them in a series expansion of Sonine polynomials and 
will consider only the leading terms. The procedure is described in Appendix IB] and only the final expressions will be 
provided here. 



B. Mass flux 



In dimensionless form, the transport coefficients associated with the mass flux, D, Dp, and D' can be written as 

nT 



D = ^^ D*, Dp - —D* D' = —D'* 

mi 7712 1^0 P'^O PVQ 



(49) 



where vq = nai2 ^0 is an effective collision frequency. The explicit forms are 

-1 



D* =[,.*- Ic 



d 
dxi 



a;i7i 



P,T 



'-^\ (l-^]D 



(50) 



-C* = xi 71 



X2 + pxi J \ 2 Zv* 



(51) 



D' 



C 
2v* P 



(52) 



Here, C* = C^^V'^o, and v* is given by 



dr(f) 



^T-(l + ai2) 



7lC2 



1/2 



{X2fj-21 +Xifii2) ■ 
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(53) 



Since j-^ = — J2 and Vxi = — Vx2, -D* must be symmetric while D* and I?'* must be antisymmetric with respect 
to the exchange 1 ^^ 2 . This can be easily verified by noting that a;i7i + 0:272 = 1- The expressions for D* , D* and 
D'* reduce to those recently obtained^^ in the tracer limit {xi -^ 0). 



C. Pressure tensor 



The shear viscosity coefficient t] can be written as 



V= — (a;i7?'7i +2;27l'72) 



where the expression of the (dimensionless) partial contribution 77* is 

72(2r22 -C*) -27iri2 



^ 7i72[C* - 2C*(tii + T22) + 4(Tiir22 - T12T21)] ' 
Here, we have introduced the (reduced) collision frequencies th and T12 given by 



(54) 



(55) 



Til = 



d-1 



A — ] (20i)"'/'(3 + 2d-3an)(l+an) 



27r('^-i)/2 

rf(d + 2)r(|) l"^' V<^i2 

+2x2^21(1 + a^2)eT0-^'^ [{d + 3){^iu02 - f^2i0i)eA{di + 62)-^'^ 

+ ^ + %'^"" m2A-^(^i + e2Y'^ + ^^|^^^r^(^i + ^.)-^/^' 



(56) 



ri2 



^-''-'^" ,^21,3/2,-1/2. ,^ . 



d(d + 2)r(|) M12 
(d + 3)(Mi2e2 - f^2i0i)9^^{ei + 62)-'^^ 

2d(d+l)-4^_, 



3 + 2d-3ai2 2/ 

-M2l6'2 ( 



^1+02)^/' 



2 ^^^"^ ^"^ ' "^' 2(d-l) "2 V^i 

A similar expression can be obtained for 772 by just making the changes 1 ^^ 2. 



-1/2 



(57) 



D. Heat flux 



The case of the heat flux is more involved since it requires to consider the second Sonine approximation. The 
transport coefficients appearing in the heat flux, D" , L, and A can be written as 



D" = - 



d + 2 



2 (mi + m2)vo 



^di + ^d*-i:^- ^' 



/il2 



/i21 



Ml2 M2I 



D* 



(58) 



d + 2 



T 



2 {1711+1112)1^0 



^^* + ^^; 



M12 



M21 



A = - 



d + 2 



nT 



2 (mi + 1112)1^0 



^Xl + ^X* 



M12 



/i21 



71 


72 


/il2 


M2I 


71 


72 


Ml2 


M21 



d: 



D' 



(59) 



(60) 
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where the coefficients D*, D* and D'* are given by Eqs. (|50l) - ((52| . respectively. The expressions of the (dimensionless) 



Sonine coefficients d*, £*, and A* are 



dl = ^{2 [2^/12^2 -yi(2;^22-3C)]Kj^2i-t^iii^22 + 2^1 + i.22)C-4C'] 



+2 f 1^ j (Fs + n) [2i^i2i^2i + 2:^22 - C*(7j^22 - 6C*)] 
{Yi + Ye) {2,^n + 21^22 -7C)\, 



-2i^l2 



dC 
dxi 



P,T 



(61) 



— {~2Y-i [2(t/i2J^21 - t^lll'22)i^22 + C{TviiV22 - 5i^l2l'21 + 2l^. 

+2Kii/i2 [21/12:^21 - 2i/iiJ/22 + 2C*(i^ii + 1^22) - C*^] 

2Y^C [2^121^21 + 1^22(2^^22 - 7C*) + 6C*'] - 2ui2CYe [2(z/ii + 



6i/iiC-7i/22C+6r')] 



1^22 



-7C1} 



(62) 



where 



K = 



1 



- {-FaC* [2^12^21 



j.22(2;^22 - 7C*) + 6C*'] + i^i2C*n [2{vii + 1^22) - 7C*] 

-n [4i^i2J^2i(j^22 - C) + 2i^2'2(5C* - 2i/n) + 2z.n(7;/22C* - 6C*') + 5C*'(6C* 
+1/12^6 [41^121^21 + 2z^ii(5C* " 2^.22) + C*(10t^22 - 23C*)]} , 



71^22)] 



A = [4(:/i2l^21 - I'lli^22) + 6C*(l^ll + J^22) - 9C*^] [l^l2J^21 " i^lli^22 + "^C {^11 + ^22) " 4C*^] 



(63) 



(64) 



In the above equations, the Y's are defined by Eqs. (jB27p - (jB29p . while the (reduced) collision frequencies vn and 
1^12 are given by Eqs. (JCITP and (jClSp . respectively. The expressions for djj ^21 ^-^id ^^2 can be obtained from Eqs. 
(pT|) - (|M|) by setting 1 <-> 2. As expected, our results for the heat flux show that D" is antisymmetric with respect 
to the change 1 «-^- 2, while L and A are symmetric. Consequently, in the case of mechanically equivalent particles 
(mi = TO2 = 771, (Ti = (72 = CT, Qfy = a), the coefficient D" vanishes. 

In the three-dimensional case {d = 3), all the above expressions for the transport coefficients reduce to those 
previously derived for hard spheres^ii^ when one takes Maxwellian distributions (|28p for the zeroth-order approxi- 
mations /j . For mechanically equivalent particles, the results obtained by Brey and Cubero^ for a d-dimensional 
monocomponent gas are also recovered. This confirms the self-consistency of the results derived here. 



E. Comparison with other theories 



Before checking the accuracy of our expressions by comparing them with computer simulations, it is instructive first 
to make some comparison with previous resultsi^i^i These results assume energy equipartition (T^ = T) and so, they 
are based on a standard CE expansion around the Maxwellian ([5^ instead of the local HCS distribution. Figures 
[H and [3] show the dependence of the reduced pressure diffusion coefficient Dp{a)/Dp{l) and the reduced thermal 
conductivity coefficient A(a)/A(l), respectively, as a function of the (common) coefficient of restitution a^ = a for 
d = 3, w = 1, xi = 0.2, and two different mass ratios fi: 11 = 0.5 (a) and ji = A (b). Here, Dp[\) and A(l) are the 
values of Dp and A for elastic collisions. We see that the deviation from the functional form for elastic collisions is 
quite important in both theories, even for moderate dissipation. It is apparent that the dependence of the transport 
coefficients on dissipation is quantitatively different in both models, especially at strong dissipation (say for instance, 
a — 0.5). This clearly shows the real quantitative effect of two different species temperatures on transport in granular 
mixtures. 



COMPARISON WITH MONTE CARLO SIMULATIONS 



A said before, the expressions derived in Sec. lIVI for the NS transport coefficients have been obtained by considering 
two different approximations. First, since the deviation of fl from its Maxwellian form ((28)) is quite small in the 
region of thermal velocities, we have used the Maxwellian distribution ([^5)1 as a trial function for /^ . Second, we have 
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FIG. 2: Plot of the reduced pressure diffusion coefficient Dp{a) / Dp(l) as a function of the (common) coefficient of restitution 
a for binary mixtures with xi = 0.2, oj = 1 in the case of a three-dimensional system (d = 3) and two values of the mass ratio 
/x: /x = 0.5 (a) and /i = 4 (b). The solid lines refer to the results derived here and the dashed lines correspond to the results 
assuming the equality of the partial temperatures. 
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FIG. 3: Plot of the reduced thermal conductivity coefRcient A(q:)/A(1) as a function of the (common) coefffcient of restitution 
a for binary mixtures with xi = 0.2, lj = 1 in the case of a three-dimensional system (d = 3) and two values of the mass ratio 
/i: /x = 0.5 (a) and /i = 4 (b). The solid lines refer to the results derived here and the dashed lines correspond to the results 
assuming the equality of the partial temperatures. 



only considered the leading terms of an expansion of the distribution /^ in Sonine polynomials. Both approximations 
allow one to offer a simplified kinetic theory for a d-dimensional granular binary mixture. To check the accuracy of the 
above predictions, in this Section we numerically solve the Boltzmann equation by means of the DSMC methodic and 
compare theory and simulation in the cases of the diffusion coefficient D (in the tracer limit) and the shear viscosity 

coefficient rj. Previous comparisons carried out for hard spheresi^ii^ when one takes into account the deviations of /^ 
from their Maxwellian forms have shown good agreement between theory and simulation, even for strong dissipation 
(say a > 0.5). Here, we expect that such good agreement is also maintained in the case of hard disks when one 



replaces /, 



(0) 



fi_M- Let us study each coefficient separately. 
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A. Tracer diffusion coefficient 

We consider the special case in which one of the components of the mixture (say, for instance, species 1) is present 
in tracer concentration {xi — > 0). In this situation, dC^* /dxi — > and so, the expression (|5D|) for the reduced diffusion 
coefficient D* becomes 

D* = ^^, (65) 



where now 



Jd-l)/2 



f ^2 



dr{d/2) \ai2 
' dV id/2) 



^2^(1 - al^), (66) 

M2i\/Aii2 + M2i7(l + "12)- (67) 



The diffusion coefficient of impurities in a granular gas undergoing homogeneous cooling state can be measured in 
simulation from the mean square displacement of the tracer particle after a time interval t'J^^ 

|(|r(i)-r(0)p> = ^. (68) 

Ot 112 

Equation (|68p is the Einstein form. This relation (written in appropriate dimensionless variables to eliminate the time 
dependence of D{t)) was used in Ref. ^M to measure the diffusion coefficient for hard spheres. More details on this 
procedure can be found in Ref. llSl . 

If the hydrodynamic description (or normal solution in the context of the CE method) applies, then the diffusion 
coefficient D{t) depends on time only through its dependence on the temperature T{t). In this case, after a transient 
regime, the reduced diffusion coefficient D* — {mim2/ p)D{t)vQ{t)/T{t) achieves a time-independent value. Here, we 
compare the steady state values of D* obtained from Monte Carlo simulations with the theoretical predictions given 
by the first Sonine approximation (|65[) . The dependence of D* on the common coefficient of restitution aij = a is 
shown in Fig. [4] in the case of hard disks for three different systems. The symbols refer to computer simulations while 
the lines correspond to the kinetic theory results given by Eq. ([65|) . Molecular dynamics (MD) results reported in 
Ref. [2y when impurities and particles of the gas are mechanically equivalent have also been included. We observe 
that MD and DSMC results ioi jjl — u = 1 are consistent among themselves in the range of values of a explored. 
This good agreement gives support to the applicability of the inelastic Boltzmann equation beyond the quasielastic 
limit. It is apparent that the agreement between the first Sonine approximation and simulation results is excellent 
when impurities and particles of the gas are mechanically equivalent and when impurities are much heavier and/or 
much larger than the particles of the gas (Brownian limit). However, some discrepancies between simulation an 
theory are found with decreasing values of the mass ratio mi/m2 and the size ratio cti/(T2- These discrepancies are 
not easily observed in Fig. 2] because of the small magnitude of D* for /i = 1/4. The above findings agree with 
those previously reported for hard spheres, ^^ where it was shown that the second Sonine approximation improves the 
qualitative predictions from the first Sonine approximation for the cases in which the gas particles are heavier and/or 
larger than impurities. The comparison carried out here for disks confirms the above expectations and shows that 
the Sonine polynomial expansion exhibits a slow convergence for sufficiently small values of the mass ratio /i and/or 
the size ratio to. 

B. Shear viscosity coefficient 

The shear viscosity 77 is perhaps the most widely studied transport coefficient in granular fluids. In the case 
of granular mixtures, this coefficient has been measuredi^ when the system is heated by the action of an external 
driving force (thermostat) that exactly compensates for cooling effects associated with dissipation of collisions. The 
corresponding shear viscosity of the mixture (which slightly differs from the one obtained in the free cooling case) has 
been determined by means of the CE method in the low-density regime^ as well as for a moderate dense mixturei^ 
The theoretical predictions compare reasonably well with the corresponding numerical solutions of the Boltzmann 
and Enskog kinetic equations. 

More recently, a new alternative method has been proposed to measure the (true) NS shear viscosity coefficient.— 
The method is based on the simple shear flow state modifled by the introduction of a deterministic non-conservative 
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D 




FIG. 4: Plot of the reduced diffusion coefficient D* as a function of the (common) coefficient of restitution a for binary mixtures 
with cj = /i in the case of a two-dimensional system {d = 2). The symbols are computer simulation results obtained from the 
mean square displacement and the lines are the theoretical results obtained in the first Sonine approximation. The DSMC 
results correspond to /x = 1/4 (•), n = A (O) and M = 1 (^)- Molecular dynamics results reported in Ref. l2a for /x = 1 (A) 
have also been included. 

force (which compensates for the coUisional cooling) along with a stochastic process. While the external force is 
introduced to allow the granular fluid to approach a Newtonian regime, the stochastic process is introduced to mimic 
the conditions appearing in the CE method to NS order. Although the method was originally devised to a single 
granular gas, its extension to multicomponent systems is straightforward. Here, we use this procedure to measure 
the shear viscosity of the mixture by means of the DSMC method. More technical details on this procedure and its 
application to dense gases can be found in Ref. [2^. 



, 1.3 




FIG. 5: Plot of the reduced shear viscosity r)*{a) = r](a)/r](l) as a function of the (common) coefficient of restitution a for 
binary mixtures constituted by particles of the same mass density (tj — fi^''^) in the case of a two-dimensional system (d = 2). 
The symbols are computer simulation results and the lines are the theoretical results obtained in the first Sonine approximation. 
The DSMC results correspond to jj, — 1 (O), ^ = 4 (•) and fi = 8 (O). We have also included DSMC results obtained in Ref. 
I29I for 11= 1 (A) from the Green-Kubo relation. 



Comparison between the first Sonine approximation and computer simulations for ri*{a) ~ iri{a)/rj{l) is shown in 
Fig. [5] for three different mixtures constituted by particles of the same mass density (i.e., /i = lo'^) in the case of a 
two-dimensional system. Here, 77(1) refers to the elastic value for the shear viscosity coefficient and we have assumed 
again a common value of the coefficient of restitution a. The symbols represent the simulation data and the lines 
correspond to the theoretical results. We have also included recent simulation results^^ for rf obtained from the 
Green-Kubo relation in the one-component case. Good agreement among the data presented here and those reported 
in Ref. [2^ for /z = 1 is observed. In addition, as happens for hard spheresji^ we see that in general the agreement 
between the first Sonine approximation and simulation is quite good. At a quantitative level, the theory slightly 
overestimates the simulation data, especially for strong dissipation and for mixtures of particles of different masses 
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and/or sizes. However, such discrepancies are quite small since for instance, they are smaller than 3% at a = 0.5 for 
/i = 8 and ui = v8. This shows again the reliability of the first Sonine approximation for the shear viscosity coefficient. 
It must be noted that this conclusion cannot in principle be extended to the transport coefficients associated with the 
heat flux since recent comparisons for a single gas^^i^i^ have shown significant discrepancies between the first Sonine 
approximation and computer simulations for high inelasticity (say a < 0.7). In this case, the agreement between 
theory and simulation can be significantly improved by the use of a modified first Sonine approximation!^ 



VI. DISCUSSION 



The main objective of this work has been to obtain the NS transport coefficients of a granular binary mixture at 
low density. In contrast to previous worksf^^ii the present study is based on a modified CE solution of the inelastic 
Boltzmann equation that takes into account non-equipartition of energy. There is no phenomenology involved as the 
equations and the transport coefficients have been derived systematically from the inelastic Boltzmann equation by 
the CE expansion around the local HCS. Since the spatial gradients are assumed to be independent of the coefficients 
of restitution, although the NS equations restrict their applicability to first order in gradients the corresponding 
transport coefficients hold a priori to arbitrary degree of inelasticity. All the calculations have been performed in an 
arbitrary number d of dimensions, previous results^^ being recovered for d = 3. 

The constitutive equations to NS order for the mass fiux, the stress tensor, and the heat fiux are given by Eqs. p3|) - 
(PS)) . respectively. The associated transport coefficients are the mutual diffusion coefficient D, the pressure diffusion 
coefficient Dp, and the thermal diffusion coefficient D' in the case of the mass flux, the shear viscosity coefficient ?/ for 
the pressure tensor, and the Dufour coefficient D", the pressure energy coefficient L, and the thermal conductivity 
A in the case of the heat fiux. These coefficients are determined from the solutions of the set of coupled linear 
integral equations (|43p - (P51) . In addition, the NS transport coefficients also depend on the reference distributions 

/j , which are not Maxwellians because they obey the integral equations ([^5]) . To solve the above integral equations 
and provide good estimates for the transport coefficients, we have considered two approximations: (i) the distributions 
/j have been replaced by their Maxwellian forms (P5|) at the temperature Tj for that species and, (ii) we have only 
considered the leading terms in a series of Sonine polynomials for the first-order distribution /j . By using both 
approximations, explicit expressions of the seven NS transport coefficients have been obtained as functions of the 
coefficients of restitution and the concentration and the ratios of mass and diameters. In dimensionless forms, the 
coefficients D, Dp, and D' are given by Eqs. ([5Dl) - ([5^ . respectively, the shear viscosity rj is given by Eqs. ([5^ and 
([55)) . while the expressions of the coefficients D" , L, and A are provided by Eqs. ((55 1) - ([M)) . 
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FIG. 6: Plot of the reduced shear viscosity coefficient rj* = ri{a)/ri{l) as a function of the (common) coefficient of restitution a 
for binary mixtures with xi = 0.2, uj — 1 and two values of the mass ratio /i: /i = 0.5 (a) and fi — 4 (b). The solid lines refer 
to spheres (d — 3) while the dashed lines correspond to disks (d = 2). 



Previous resultsS' derived from the CE method have typically introduced additional assumptions for convenience 
that are not internally consistent with constructing a solution to the Boltzmann equation. Thus, in most of the cases 
the reference state /^ has been chosen to be a Maxwellian at the same temperature [see Eq. ([5^ ]. This assumption 
is presumed to give accurate results at weak dissipation where energy equipartition can be still considered as a 
good approximation. However, as shown in Fig. [U the temperature ratio T1/T2 clearly differs from 1 as dissipation 
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increases. Here, we have replaced /^ -^ fi^^i so that, the influence of the fourth-cumulants Ci of f> has been 
ignored.!^ Comparison between the expressions derived in this paper by taking the approximation (|28p with those 
obtained by assuming energy equipartition shows important discrepancies as the coefficient of restitution decreases 
[see Figs. [2] and [3]. Moreover, as an added value of our theory, the use of the Maxwellian approximation ([28)1 for fl ' 
allows one to provide simple and explicit expressions for all the transport coefficients in terms of the parameters of 
the mixture. This contrasts with the relatively recent study for hard spheres^^ where the constitutive relations for 
the fluxes were not explicitly displayed. 

As a complementary route and to check the reliability of our theory, the analytical results derived for the diffusion 
coefficient D and the shear viscosity rj in the first Sonine approximation have been compared with those obtained 
from numerical solutions of the Boltzmann equation by means of the DSMC method for a two-dimensional system. 
For the sake of simplicity, all the simulations have considered a common coefficient of restitution a = aij. As 
expected, theory and simulation clearly show that the influence of dissipation on mass and momentum transport is 
quite important since there is a relevant dependence of the diffusion D and viscosity 77 coefficients on a. With respect 
to the accuracy of the theoretical predictions, we see that in general the CE results in the first Sonine approximation 
exhibit a good agreement with the simulation data. Exceptions to this agreement are extreme mass or size ratios 
and strong dissipation. These discrepancies are basically due to the use of the first Sonine approximation and can be 
partially mitigated by considering the second and third Sonine approximationai^ or the use of a modified first Sonine 
approximationj2E 

As said in the Introduction, the results obtained in this paper are of great practical interest since most of the 
experiments and simulations are performed in two dimensions. On the other hand, apart from this practical interest, 
the knowledge of the NS transport coefficients of a d-dimensional mixture allows one to investigate the influence of 
dimensionality on the transport properties of the system. To illustrate this effect, in Fig. [S] we plot the reduced shear 
viscosity 77* = ry(a)/?7(l) versus the coefficient of restitution a for to = 1, xi = 0.2, and two different mass ratios /z: 
fi = 0.5 (a) and IJ, — ^ (b). We have considered the physical cases of hard spheres (solid lines) and hard disks (dashed 
lines). Although the qualitative dependence of rf on a is quite similar in both systems, we observe that the influence 
of dissipation on momentum transport is stronger for rf = 3 than for rf = 2. This trend is also observed in general in 
the remaining transport coefficients. 

One of the main limitations of our theory is its restriction to dilute gases. In this situation, the coUisional transfer 
contributions to the fluxes are neglected and only their kinetic contributions are considered. Possible extension of the 
present kinetic theory to higher densities can be done in the context of the revised Enskog theory. Preliminary results- 
have been focused on the uniform shear fiow state to get directly the shear viscosity coefficient. The extension of this 
studji^l to states with gradients of concentration, pressure, and temperature is somewhat intricate due to subtleties 
associated with the spatial dependence of the pair correlations functions considered in the revised Enskog theory. On 
the other hand, it must be remarked that many of the collision integrals appearing in the Enskog description are the 
same as those appearing in the Boltzmann limit so that one can take advantage of the results reported in this paper. 
We plan to extend the results derived for moderately dense mixtures of smooth elastic hard spheres^ to inelastic 
collisions in the near future. 
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APPENDIX A: CHAPMAN-ENSKOG METHOD 



The velocity distribution function /-[ obeys the equation 



d 



(9i")+A)/P+A^i/i^^ = - a(^^+v.V + g.|- /f', (Al) 



where the linear operators £1 and M.i are defined by Eqs. (|47p and (|48p . respectively. A similar equation can be 

3ft 

d'^^^^xi = 0, (A2) 



obtained for /2 by interchanging 1 <-» 2. The action of the time derivatives d\ on the hydrodynamic fields is 



d[% 



d+2 



pVu, 



17 
(A3) 



l(l)r 



2T. 



D-'T^--V.u, 



(A4) 



A^'^u = -p-iVp + g, 



(A5) 



i(i) 



9f^) 



:(0) 



where D'j.'^' — dY' + u • V and use has been made of the results j^"^ — q'"-* — C,'^^^ = 0. The last equality follows 
from the fact that the cooling rate is a scalar, and corrections to first order in the gradients can arise only from the 
divergence of a vector field. However, as is demonstrated below, there is no contribution to the distribution function 
proportional to this divergence. We note that this is special to the low density Boltzmann equation and such terms 
do occur at higher densitiesi^ Use of Eqs. (|A2p - (|A5p yields 



a« + v.v + g.A)/<"' = 



l/'X.^-^- 






Vlnp 



f\ 



(0) 



1 d 
29V 



V/{ 



(0) 



V- VlnT 






1 



/kiV ■ -^.ff^ I VfeM^ 



(A6) 



Note that the external field does not appear in the right-hand side of Eq. (jA6|) . This is due to the particular form of 
the gravitational force. Using Eq. (|A6[) . Eq. (jAip can be written as 



^.(0) 



p(i) 



p(i) 



£1 fl^' +MifY' - Ai • Vxi + Bi • Vp + Ci • VT + Di : Vu, 



(A7) 



where 



A.(V) 



dXl Jprj. 



(A8) 



B.(V) 






(A9) 



C.(V) = - 



1 d 






V, 



(AlO) 



D.(V) 






iiv.A/(o). 



(All) 



In Eqs. (|A8|) - (|Alip it is understood that i ~ 1,2 and 1 is the unit tensor in d dimensions. Note that the trace of D^ 
vanishes, confirming that the distribution function does not have contribution from the divergence of the flow flcld. 
The solutions to Eqs. (|A7p are of the form 



/f ^ = A, ■ Vxi + B, ■\/p + C^-\IT + I?,,fcfVfcWf 



(A12) 



The coefflcients A^, Bi, Ci, and X>i are functions of the peculiar velocity V and the hydrodynamic flelds. The cooling 
rate depends on space through its dependence on xi, p, and T. The time derivative d^ acting on these quantities 
can be evaluated by the replacement d^ — > — (■(") {Tdr ' "^^ '^ t„ „^^u;^„ +u„v.„ „,.„ „„^+,.;k,,+;^v,o r,,^^ ay ) 
on the temperature and pressure gradients given by 



■pdp). In addition, there are contributions from d^. acting 



a^°VT 



-V (re'" 



-(^("'VT-TVC 



(0) 



\/T -T 



\dxj, 



Vxi 



(^(0) 



Vp 



(A13) 
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^2C^"^Vp-p 



dxi 



P-T 



HO) 
Vxi - ^ — VT 
2T 
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(AM) 



The corresponding integral equations for the unknowns jA-t, Bt, Ct, and "Di are identified as the coefficients of the 
independent gradients in (|A12p . This leads to Eqs. (|i5)) " ((i^ . 

APPENDIX B: LEADING SONINE APPROXIMATIONS 

In this Appendix, we get the explicit expressions of the mass, momentum, and heat fluxes in the first Sonine 
approximation and neglecting the non-Gaussian corrections to the reference distributions /^ (i-G-, the cumulants 
Ci = 0). The procedure to get the leading order contributions in the Sonine polynomial expansion to the transport 
coeflficients is quite similar to the one previously used in the three-dimensional case.— Only some partial results will 
be presented here. 



1. Leading Sonine approximation to Mass Flux 



In the case of the mass flux, the leading Sonine approximations (lowest degree polynomial) of the quantities ^t, 
£5j, K^i are 



A(V)^-/,Mv!!^i!^i?, A,iV) ^ f,^,,Y!!^l^D 



pniTi 



pn2T2 



(Bl) 



Bi(V) ^ -/i,mV— ^i^p, 02 (V) ^ f^^^tY^-Dp 



(B2) 



Ci(V) -. -f,^^,V-J^D', C2(V) ^ /2,MV-ii-i?', 
Inili in2i2 



(B3) 



where fi^M are the Maxwellian distributions ([^5]) . Multiplication of Eqs. (|^5|) - P5)) by toiV and integrating over the 
velocity yields 



C'-"HTdT+pdp) + iy 



mim2n 



D 



^.-'\,~^^-£)j'^-^''^ 



(B4) 



C''^{TdT+pdp)-2C'^''^ + 



(0) 



'-D, 



niTi f _ minT \ _ PC^jji 
pTi 



(B5) 



Here, v is the collision frequency defined by 



2T ^ 



l-JdY,m,V, ■ [£i(/i,mVi) - <57Mi(/2,A/V2)] 
^^JdV.rmV, ■ (Ji2[vi|/i,mVi,/("^] - <57Ji2[vi|/r, /2,mV2 



(B6) 



(B7) 



where S = X\lx2- The evaluation of the collision integral (|B7p is made in Appendix [Cl The self-collision terms of Li 
arising from Jn do not occur in Eq. (|B7[) since they conserve momentum for species 1. From dimensional analysis, 
D ~ T"'^/^, Dp ~ pT~^^'^, and D' ~ pT~^l'^ so the temperature and pressure derivatives can be performed in Eqs. 
(jB4p - (jB6l) . After performing them, one gets the expressions ([501) . I|5ip . and (|5^ for the (reduced) coefficients Z?*, 
Z?* and £>'*, respectively. 
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2. Leading Sonine approximation to Pressure Tensor 

In the case of the pressure tensor, the leading Sonine approximation for the function Vi ki is 



where 



and 



V^Mi^) -> -/.,M(V)|i?.,M(V), i = l,2 



R^Mi^)^^^ [VkVi--V^6ki ), 



1 T 



'' {d~l)id + 2)n,Tl 

The shear viscosity 77 in this approximation can be written as 



y"dvi?,,fe,(V)2?,,fe,(V). 






(B8) 



(B9) 



(BIO) 



(Bll) 



where ry* = I'o'^i- The integral equations for the (reduced) coefficients ry* are decoupled from the remaining transport 
coefficients. The two coefficients 77* are obtained by multiplying Eqs. PS|) with Ri^M and integrating over the velocity 
to get the coupled set of equations 



ri\ _ 7i 



Til - 2C* '''12 _ 

T21 T22 -\C )\ V2 J V 72^^ 



The (reduced) collision frequencies Tij are given in terms of the linear collision operator by 



1 



(d-l)(d + 2)n,r, 



-2 — / dviRikid {fi,MRi,ke) ■ 
-i ^0 J 



(B12) 



(B13) 



1 



{d-l){d + 2)n,T^ 



^ — / dviRi,keMi Uj.mRjm) , i 7^ J- 



(B14) 



The evaluation of these collision integrals is also given in Appendix [Cl The solution of Eq. (|B12p is elementary and 
yields Eq. ([55)) . 



3. Leading Sonine approximation to Heat Flux 



The heat flux requires going up to the second Sonine approximation. In this case, the quantities ^, Bi, Ci are 
taken to be 



Ai{\)^h. 



M 



pniTi 



-DY + d'(SiiY) 



^2(V) ^ /2, 



M 



mim2n 



pn2T2 



-DV + 4'S2(V) 



(B15) 



Bi(V) ^ /i 



M 



pniTi 



■DpV + iiSiiY) 



B2(V) ^ /2. 



M 



pn2T-^ 



■DpV + e2S2{V) 



(B16) 



Ci{V) ^ /i. 



M 



TniTi 



D'V + XiSiiV) 



, C2(V)^/2, 



M 



Tn2T2 



i:)'V + A2S2(V) 



where 



1 X.2 d + 2^ 
-niiV 

2 2 



S.(V) = ( -m.1/2 - -^-^Ti \ V. 



(B17) 



(B18) 
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In these equations, it is understood that D, Dp and D' are given by Eqs. ([50| . (I5ip . and ((52| . respectively. The 



cocfBcients d", ^^ and A,; are defined as 



t 



^_ I d(d + 2)n. 



»/v^ 



^ /dvS,(V). S, 



(B19) 



These coefficients can be determined by multiplying Eqs. (j43|) -([45 |l (and their counterparts for the species 2) by Si(V) 
and integrating over the velocity. The final expressions can be obtained by taking into account that d![ ^ T^'^l'^ ^ 
li ~ T-^^^/p, and Ai ~ T'^/^ and the results 



/'dvmiSi(V)-Ai 



-7i 



nil \dxi 



P,T 



(B20) 



dvmiSi(V)-Bi = 0, 



(B21) 



/"dvTOiSi(V)-Ci = - 



d{d + 2) n^Tl 
2 miT' 

By using matrix notation, the coupled set of six equations for the quantities 

{did;,iii;,\ix*^} 

can be written as 



(B22) 



(B23) 



(B24) 



Here, d* = Tv^d'l , i* = pTv^li, A* = T^i^oXi, Xa-' is the column matrix defined by the set (|B23|) and Ag-cr' is the 
square matrix 



/ ^11 - ic 



A = 



i^l2 







V 



.(if.) -(f) \ 

I'll - oC 



1^12 



5 /"* 
!^21 1^22 - 2C 



CV2 
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i^21 





-C 

i^l2 
i^22 - C* / 



The column matrix Y is 



wher( 



,36 



Yi 



icoi2 - a 



1 /971 



D* 



xi"il 7i \dxijpj.' 



i"2 
X3 



Fo 



D* 



X27J 



2 (W21 - C ) - -2 



7I V^a^iyp 



(B25) 



(B26) 



(B27) 



D* D* 
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1 D'* 1 D'* 

5-5 = + 2 (^12 - C) , ^6 = 2 (^21 - C) ■ 

71 a^i7i 72 a;272 



(B28) 
(B29) 
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Here, we have introduced the (reduced) coUision frequencies 

2 



d{d + 2) 



Tfl ■ f 

— — = — / dviSi • Ci [fiM^i) , 
niTfiyo J 



d{d + 2) 



Tfl ' f 

— -|— / dViS, • M^ (/j,A/Sj) , l^j, 



(B30) 



(B31) 



Wl2 



mi 



d{d + 2) niT^iyQ 



/dviSi • /:i(/i,AfVi) - ^7 /dviSi • A1i(/2,mV2) 



(B32) 



t^21 = 



m.2 



dviSa • C2if2.M'Vi) - T- / dviSa • 7W2(/i,mV2) 



(B33) 



d{d + 2) n2T^iyo 
The expressions of the colhsion integrals (|B30[) . (|B3ip . and (|B32p are given in Appendix [Cl The solution to Eq. (|B24p 



X,^{A-')^^,Y„,. (B34) 

From this relation one gets the expressions (pTj) . ([5^ . and ([55| for the coefficients d*, £* and A*, respectively. 



APPENDIX C: COLLISION INTEGRALS 

In this Appendix we compute the different collision integrals appearing in the expressions of the transport coeffi- 
cients. To simplify all the integrals, we use the property 

dvM^l)J^3 [vil/.,/,] = Tfi' I dvi / dv2/,(Vi)/,(V2) t da Q{a ■ gi2)(a • gi2) \h{V[) - h{YA , (Cl) 



with 



Vi == Vi - iiji{l + aij){a- ■ gi2)o^ 



This result applies for both i = j and i ^ j. 

Let us start with the collision frequency v defined by Eq. (jB7p . Use of the identity (IC2[) in (jB7|) gives 



V = ^^53^i2"V2l(l+ai2) / dVi / dV2<?12f/l,M(Vi)/2(°)(V2)(Vi-gi2) 



dniTi 

-'57/l"^(Vi)/2,M(V2)(V2-gl2)' 



where use has been made of the result 



/ da e{a ■ gi2) {a ■ gi2)'=? = Bfc+iffi2 ^gi2, 



whith^^ 



Bk= da e{a ■ gu) (a ■ g^f = tt 



f(0) 



k _ „(d-l)/2^ I 2 ) 



r(^)' 



Substitution of the Maxwellian approximation (|28p for fl gives 
2 7r(''-i)/2 



dr(^) 



■^^M^ + ai2)7r-'' iei02f^ J dci y dc2 2/e-(^^^?+''^^^) [a;27r'(ci • y) - ^172-^02 • y)] , 



(C2) 



(C3) 



(C4) 



(C5) 



(C6) 



22 

where c^ = Vi/uo and y = Ci — C2. The integral can be performed by the change of variables {01,02} ^ {y, z}, where 
' 02 C2 and 

2 7r('^-i)/2 



z = 9iCi + 92C2 and the Jacobian is {61 + O2) ■ With this change, Eq. (|C6P becomes 



v = — 



dT(^' 



-M^ + «i2)^-' (01^2)^'+'^/' {Oi + 92)-^'+"^ (X2A*21 + .TiA*i2) / dy / dzy3e-(-^'+^-'), (C7) 



where a = 6'i6'2 (6*1 + (^2) and b = {9i + 62)^ ■ The integral (jC7p can be easily computed and one directly gets the 
result (j53p given in the text for the reduced collision frequency v* — v /vq. 

The collision frequencies r^ defined by Eqs. (|B30P and (jB31|l involve collision integrals of the form 

ydviViVij,,[/„/,] -^-^y dviy dv2/.(vi)/,(V2)ydae(^-gi2)(^-gi2) [v/v/- vivi], (cs) 



where the identity (jCip has been used. The scattering rule (jC2p gives 

V'/V" - ViVi = -^ji(l + aij){o- ■ gi2) [GijCT + a-Gij + /Xji(gi2CT + CTgi2) - ^J^.ji{l + aij){d- ■ gi2)?CT] , (C9) 



where Gy = fJ-ij^i +/^jjV2. Substitution of Eq. (|C9|1 into Eq. (jCSP allows the angular integral to be performed with 
the result 

y da e{a ■ gi2)(S^ • gi2) [V'/V/ - ViVi] = -B3m,/x,,(l + ay) [5i2(G,,gi2 + gi2G,y) 

+3^(l + y-«^j)5i2gi2gi2-^(l + ay)/l . (CIO) 

Using (|C10|) the integrals defining r.y can be calculated by the same mathematical steps as those made before for v. 
After a lengthy calculation, one gets 

/ d-viRi^keJi2[.fi , f2,MR2M] = --r-7r-7rfY"^i'^2»^i"-2A*2i(l + ai2)cri2"^Wo (^1^2)" 

X {2(d + 3)(d - I){^li2e2 - fl2l0l)9^^ (01 + 02)^'/' 

+3(d - i)m2i (1 + y - «12) ^2"' C^i + ^2)'/' 

- [2d(d + 1) - 4] 9-' (01 + 02)"'^'} , (Cll) 

/ dvii?i,MJi2[/l,Jl/i?l,fc£,/2 ] == --y-r^ "117712711712^^21 (1 +ai2)o-i2'^t'0 (0102)" 

X {2(d + 3)(d - 1)(a*i202 - /i2i0i)0r' (01 + 02)"'/' 



+3{d - l)fi2i ( 1 + y - ai2 ) 0r' (01 + 02)'/' 



+ [2d{d + 1) - 4] 9-^ i9i + 02)"'/'} , (C12) 

Jdv^RiM {jii[fi°\fiMRiM] + Jii[fiMRiMJi°'']} = - ""^ri) "^1^1(1 + "ii)'^i"'(^iMi)'/' 



.6{d~l}r 2d 



d V ■ 3 ""j ■ (^^^) 

The corresponding expressions for Tij can be easily inferred from Eqs. (JC11|) - (|C13[) . 

The collision frequencies Vij and cuij that determine the heat flux are defined by Eqs. (|B30p . (jB3ip . and (|B32[) . 
respectively. To evaluate these collision integrals, one needs the partial results 

S,(Vi ) - S,(Vi) = ^(1 + a.,,)f,,,ia ■ gi2) { [(1 - 4)4(a • gi2)' - Gl ~ f,%9l2 
-2^jj(gi2 • G,j) + 2(1 + a,j)^iji{a ■ gi2)(CT • G^) + (d + 2)—^ 

TTli 

- [(1 - a,j)^ij,{a- ■ gi2) + 2(S • Gjj)] G.y 

-liji [(1 - a^j)fiji{a- ■ gia) + 2(ct • G^j)] gi2} , (C14) 
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JdaQ{a ■ gi2)(S • gi2) [S,(V;') - S,(Vi) 
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ff?2 
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''J 
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(d + 5)ay — d — 1 3 



rrii 
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2512 (gl2 ■ Gy) -/^tji- 



d + 3 



-5i2 



(C15) 



The integrals w^- and Vij can be explicitly evaluated by using (|C15|) and the same mathematical steps as before. After 
a lengthy algebra, one gets 



T,(d-l)/2 2 



W12 



r(|) d^/2\(Jl2 
7r('^-i)/2 2 

^ r(|) d(d + 2) 
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72(,Pi -|-(72j 
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In the above equations we have introduced the quantities^ 

A ^ {d + 2)(2/3i2 + 62) + M2i(0i + 02) {(d + 2)(1 - ai2) - [(11 + d)ai2 - 5d - 7]/3i20r'} 

+3(rf + 3)/3220-i + 2/^2^ ( 2a22 _ ^^aia + rf + 1 ) 0r'(0i + ^2)' - (rf + 2)02^7' (^i + ^2), 



B = (d + 2)(2/3i2 - 0i) + ^i2l{0l + 02) {(d + 2)(1 - ai2) + [(11 + d)ai2 -5d- 7]f3i202^} 
-3(d + 3)/3220-i - 2/^2^ ( 20^2 - ^^ai2 + d + 1) 0^'(0i + 62)^ + {d + 2)(0i + ^2), 



(C18) 



(C19) 



(C20) 



E = 2nl^e^^{ei + 62)^ (2al2 - ^-y^"i2 +d+l][{d + 2)01 + {d + 5)^2] 

-^i2l{0l + 62) {/3i20r'[(rf + 2)^1 + (d + 5)02][(11 + d)ai2 -5d-7] 

-6l26lj-i[20 + d(15 - 7ai2) + ^^(1 - ^12) - 28ai2] -{d + 2)^(1 - aia)} 

+3(rf + 3)/3?20r^[(c; + 2)6ii + (d + 5)6*2] + 2/3i20rM(c^ + 2)^6*1 + (24 + lid + ^2)612] 

+ {d + 2)02e^^[{d + 8)01 + {d + 3)02] -{d + 2)(0i + 02)0r'02[(d + 2)01 + {d + 3)02], 



(C21) 
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F = 2nl^9~'{9i + 02^ ( 2al2 - ^-^a^ +d+l][{d + 5)0i + {d + 2)9, 



2 



2 

-A*2i(^i + O2) {I3i202^[{d + 5)01 + {d + 2)02][(11 + d)ai2 -5d~7] 

+6ii6i2"i[20 + rf(15 - 7ai2) + d^{l- a^) - 28ai2] + {d + 2)2(1 - 0,12)} 

+3(d + 3)/3?2^^^ [(o' + 5)6*1 + (d + 2)612] - 2/3i26i2"^ [(24 +lld + d^)ei + {d + 2fe2] 

+{d + 2)9ie:^^[{d + 3)91 + {d + 8)02] -{d + 2)(9i + 92)9:2^[(d + 3)6*1 + (d + 2)6*2]. 

(C22) 

Here, /3i2 = ^126*2— M2i6'i- From Eqs. (|C16|) - (|C22p . one easily gets the expressions for UJ21, 1^22 and 1^21 by interchanging 
1 ^ 2. 

In the case of a three-dimensional system (d = 3) , all the above results reduce to those previously obtained for hard 
spheres when one takes Maxwellian distributions for the reference homogeneous cooling state.— i" 
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